#ifndef NCS_H
#define NCS_H

#include <cassert>
#include <algorithm>
#include <vector>

#include "Curves/Interfaces/icurve.h"

using std::vector;

template< class IPoint >
class NCS : public ICurve< IPoint >
{
private:
    //vector<QPointF> control_points;
    vector< double > h, l, dx, dy,
                     p, q, ux, uy,
                     Mx, My, t;
    // void computeHelper( const vector<double> );
    void init_t( );

public:
    NCS( );
    NCS( const IPoint &point );
    NCS( const vector< IPoint > &points );
//    NCS( const NCS &ncs );
    ~NCS( );

    vector< IPoint > ComputePoints( int n );
    IPoint Eval( double t ) const;
};

template< class IPoint >
NCS< IPoint >::NCS( )
{
    this -> degree = 3;
}

template< class IPoint >
NCS< IPoint >::NCS( const IPoint &point )
{
    this -> control_points.push_back( point );
    this -> degree = 3;
}

template< class IPoint >
NCS< IPoint >::NCS( const vector< IPoint > &points )
{
    this -> control_points = vector< IPoint >( points );
    this -> degree = 3;
}

template< class IPoint >
NCS< IPoint >::~NCS( )
{ }

template< class IPoint >
vector< IPoint > NCS< IPoint >::ComputePoints( int m )
{
    if ( this -> computed_points_up_to_date )
        return this -> computed_points;

    const vector < IPoint > & f = this -> control_points;
    int n                       = this -> control_points.size( ) - 1;

    init_t( );

    p  = vector< double >( n+1, 0 );
    q  = vector< double >( n+1, 0 );
    l  = vector< double >( n+1, 0 );
    dx = vector< double >( n+1, 0 );
    dy = vector< double >( n+1, 0 );
    ux = vector< double >( n+1, 0 );
    uy = vector< double >( n+1, 0 );
    Mx = vector< double >( n+1, 0 );
    My = vector< double >( n+1, 0 );
    h  = vector< double >( n+1, 0 );

    for ( int k = 1; k <= n; ++k )
    {
        h[k] = t[k] - t[k-1];
        l[k] = h[k] / (h[k] + h[k+1]);
        dx[k] = 6.0 / (t[k + 1] - t[k - 1])
                * (   (f[k + 1].x( ) - f[k].x( )) / h[k + 1]
                    - (f[k].x( ) - f[k - 1].x( )) / h[k]     );
        dy[k] = 6.0 / (t[k + 1] - t[k - 1])
                * (   (f[k + 1].y( ) - f[k].y( )) / h[k + 1]
                    - (f[k].y( ) - f[k - 1].y( )) / h[k]     );
    }

    for ( int k = 1; k <= n-1; ++k )
    {
        p[k] =  l[k] * q[k - 1] + 2;
        q[k] = (l[k] - 1) / p[k];
        ux[k] = (dx[k] - l[k] * ux[k - 1]) / p[k];
        uy[k] = (dy[k] - l[k] * uy[k - 1]) / p[k];
    }

    Mx[n-1] = ux[n-1];
    My[n-1] = uy[n-1];

    for ( int k = n-2; k >= 1; --k )
    {
        Mx[k] = ux[k] + q[k] * Mx[k+1];
        My[k] = uy[k] + q[k] * My[k+1];
    }

    this -> computed_points.resize( n );
    for ( int k = 0; k < m; ++k )
    {
        this -> computed_points[k] = Eval( t[k] );
    }

    this -> computed_points_up_to_date = true;
    return this -> computed_points;
}

template< class IPoint >
IPoint NCS< IPoint >::Eval( double T ) const
{
    const vector< IPoint > &f = this -> control_points;
    int n                     = this -> control_points.size( ) - 1;

    int k = 1;
    // find section t[k-1] < T < t[k]
    for ( ; k < n + 1; k++ )
        if ( t[k - 1] <= T && T <= t[k] )
            break;

    double Sx = (    1/6 * Mx[k - 1] * pow( t[k] - T    , 3 )
                   + 1/6 * Mx[k]     * pow( T - t[k - 1], 3 )
                     + ( f[k - 1].x( ) - 1/6 * Mx[k - 1] * h[k] * h[k] ) * (t[k] - T)
                     - ( f[k].x( )     - 1/6 * Mx[k]     * h[k] * h[k] ) * (t[k - 1] - T)
                ) / h[k];
    double Sy = (    1/6 * My[k - 1] * pow( t[k] - T    , 3 )
                   + 1/6 * My[k]     * pow( T - t[k - 1], 3 )
                     + ( f[k - 1].y( ) - 1/6 * My[k - 1] * h[k] * h[k] ) * (t[k] - T)
                     - ( f[k].y( )     - 1/6 * My[k]     * h[k] * h[k] ) * (t[k - 1] - T)
                ) / h[k];

    return IPoint( Sx, Sy );
}

template< class IPoint >
void NCS< IPoint >::init_t( )
{
    int n = this -> control_points.size( ) - 1;

    if ( t.size( ) < static_cast< unsigned int >( n + 1 ) )
    {
        for ( double i = 0.0; i <= n; i += 1.0 )
            t.push_back( i + 1 );
    }
}

#endif // NCS_H
